QTL analysis and candidate gene prediction for seed density per silique by QTL-seq and RNA-seq in spring Brassica napus L.

Seed density per silique (SD) is an important agricultural trait and plays an important role in the yield performance of Brassica napus L. (B. napus). In this study, a genetic linkage map was constructed using a double haploid (DH) population with 213 lines derived from a cross between a low SD line No. 935 and a high SD line No. 3641, and a total of 1,098,259 SNP (single-nucleotide polymorphisms) markers and 2,102 bins were mapped to 19 linkage groups. Twenty-eight QTLs for SD were detected on chromosomes A02, A04, A05, A09, C02, C03, C06, and C09 of B. napus, of which eight QTLs were on chromosome A09 and explained 5.89%-13.24% of the phenotypic variation. Furthermore, a consistent QTL for SD on chromosome A09, cqSD-A9a, was identified in four environments by QTL meta-analysis, explaining 10.68% of the phenotypic variation. In addition, four pairs of epistatic interactions were detected in the DH population via QTL epistasis analysis, indicating that SD is controlled not only by additive effects but also by epistatic effects that play an important role in spring B. napus., but with little environmental effect. Moreover, 18 closely linked SSR markers for cqSD-A9a were developed, as a result, it was mapped to a 1.86Mb (7.80–9.66 Mb) region on chromosome A09. A total of 13 differentially expressed genes (DEGs) were screened in the candidate interval by RNA-seq analysis, which were differentially expressed in buds, leaves and siliques both between and siliques both between two parents and two pools of extremely high-SD and low-SD lines in the DH population. Three of 13 DEGs were possible candidate genes that might control SD: BnaA09g14070D, which encodes a callose synthase that plays an important role in development and stress responses; BnaA09g14800D, a plant synaptic protein that encodes a membrane component; and BnaA09g18250D, which is responsible for DNA binding, transcriptional regulation, and sequence-specific DNA binding and is involved in the response to growth hormone stimulation. Overall, these results lay a foundation for fine mapping and gene cloning for SD in B. napus.


Introduction
Oilseed rape (Brassica napus L.) (AACC, 2n = 38) is an allotetraploid species formed by natural hybridization derived from genome doubling between Brassica rapa (AA, 2n = 20) and Brassica oleracea (CC, 2n = 18) and originated in Europe [1]. B. napus produces edible oil, and its cake is rich in protein and can be used as animal feed. This species is also a raw material for industrial lubricants and biofuels. B napus is not only one of the important oil crops in the world but also one of the four important oil crop species in China (which are rape, soybean, peanut, and sesame), and rapeseed oil accounts for approximately 40% of the edible vegetable oil in China [2,3]. Therefore, improving rapeseed yield remains an important research topic.
Seed yield is a complex trait that is influenced by various factors. Thousand-seed weight (TSW), number of seeds per silique (SPS), effective siliques per plant (ESP), and silique length (SL) are closely associated with seed yield improvement. The seed density per silique (SD) refers to the number of seeds per unit silique length, that is, the ratio of the number of seeds per silique to the silique length. Studies have shown that SL and is extremely significantly positively correlated with TSW, and SD is significantly negatively correlated with TSW and SL [4, detected on chromosomes A3, A7, A9 and C3; among these QTLs, two major loci, cqSD-A9-c and cqSD-A9-d, were mapped in 7 environments [26].
There are few research reports on SD, and only a few related QTLs have been found; however, several candidate functional genes have been identified. Functional analysis and verification of these candidate genes and the cloning of key genes have not yet been reported. In addition, previous studies on SD have been limited to winter oilseed rape, and there are almost no related studies on spring oilseed rape. In this study, a 213 DH population was constructed from two spring oilseed rape lines with significant differences in SD. They were resequenced and subjected to multiyear and multipoint phenotyping experiments to identify major QTLs associated with SD. Then, the candidate genes associated with the major QTL were identified by transcriptome analysis. This will provide guidance for gene aggregation in the development of high-yield resources in the future.

Materials
Two spring oilseed rape varieties with excellent agronomic characteristics were selected. No. 935, which has a low SD, was crossed with No. 3641, which has a high SD, to produce an F 1 generation. Two thirteen DH lines were subsequently obtained by microspore culture. According to individual phenotype differences, we selected 30 DH lines with extremely high SD and 30 DH lines with extremely low SD to construct two gene pools (high SD pool and low SD pool) for QTL-Seq analysis by the BSA method. Buds, leaves, siliques, and two bud pools constructed from 20 DH lines with extremely high SD and 20 DH lines with extremely low SD (high SD pool and low SD pool) were used to construct transcriptome libraries.

Linkage map construction
The DNA of each sample was constructed for paired-end (PE) library construction after passing the quality inspection. Then, PE150 sequencing was performed on an Illumina Hi Seq sequencer. The clean data obtained through data quality control were compared to the reference genome Darmor (v4.1) (downloaded from: http://www.genoscope.cns.fr/brassicanapus/ data/) for mutation detection. (software: Cutadpat and Trimmomatic). The source of each allele is determined according to the parental genotype (parental P1 is marked as A, parental P2 is marked as B). The genetic linkage map was ultimately obtained using the hidden Markov model (HMM) algorithm [27] and the Kosambi mapping function to calculate the genetic map distance between markers. The construction of the genetic linkage map was performed by staff at Wuhan Kinosec Technology Co., Ltd.

Field trials and phenotypic investigations
The DH population, parents, and F 1 were planted at two experimental sites, Huzhu (Qinghai Province, altitude:2,557m, 101˚57 0 N, 36˚49 0 E) and Xining (Qinghai Province, altitude 2,225 m, 101˚49 0 N, 36˚34 0 E), in mid-April 2019 and early-April 2020, respectively, which were designated 2019HZ, 2019XN, 2020HZ and 2020XN. A randomized complete block design was conducted with three replications. Each plot was planted in three rows, with 10 plants in each row and a distance of 15 cm between plants within each row and 30 cm between rows. At the mature stage, 6-10 open-pollinated plants from each plot were randomly selected to test agronomic traits, and heritability. For agronomic traits, the main tests included effective silique number of plant (ESP), thousand seed weight (TSW), silique length (SL), silique length (SL), seed density per silique (SD) and yield per plot (YP). The data of these traits were analyed by using SPSS 18.0 software [28]. (H 2 : broadsense heritability, V G : genetic variance, V P : phenotypic variance, V E : environmental variance,  V A : additive variance, V D : dominant variance.) QTL analysis of SD QTL mapping was performed by the composite interval mapping (CIM) method using Win QTL Cart 2.5 according to the descriptions of Doerge et al. [29]. The 1000 permutation test method was used to determine the threshold of the LOD value for the phenotype data of each environment with p = 0.05. 10 cM window size and 1 cM walking speed are selected. QTLs detected in different environments were considered to be consensus QTLs if they overlapped, and QTL naming followed the methods of McCouch et al. [30].

QTL meta-analysis of SD
QTLs with overlapping regions detected in different environments or populations were integrated using Bio Mercator 4.2 software [31]. A consensus QTL must contain at least two component QTLs. Compared with component QTLs, the confidence interval of each integrated consensus QTL was significantly reduced, and the phenotypic variation rate was the average of the component QTLs. The QTL integrated by the model with the smallest Akaike information criterion (AIC) value is considered to be the true QTL [32].

Epistasis analysis of SD
Epistasis analysis was performed using QTL Network 2.0 software [33]. Additionally, 500-replicate permutation analysis of phenotypic data was performed in triplicate for each environment to estimate the threshold for LOD values via 2D genome scanning. The LOD threshold candidate interval and the QTL effect were determined with P = 0.05. The full QTL model was used to construct an epistatic QTL map of epistatic effects and epistatic effect× environment interactions.

QTL sequencing (QTL-seq) analysis for SD
Genomic DNA of green leaves from each phenotype was pooled in equal amounts from 60 DH lines (30 extremely high SD and 30 extremely low SD lines) to construct a high SD (HSD) pool and a low SD (LSD) pool, respectively. No. 935 and No. 3641 were used as parental pools. The libraries were sequenced using the Illumina HiSeq™ PE150 platform, and DNA samples were randomly disrupted to obtain fragments of approximately 200-500 bp with a Covaris crusher. Clean reads with a length greater than 50 bp are obtained through Cutadapt software (version 1.13) and Trimmomatic software (version 0.36) s/to remove adapter sequences and low-quality bases. Sequencing data were aligned to a reference genome using the MEM algorithm of BWA software (version 0.7.15-r1140) (Darmor v4.1 2, available at http://www. genoscope.cns.fr/brassicanapus/data/). SNPs and Indels are detected and screened after reads are compared with the reference genome. The Δ(SNP index) was calculated and used to detect the allele frequency between the HSD pool and LSD pool, and then the interval with a significant allele frequency difference between the two pools was considered the candidate interval for the QTL of SD. The stronger the correlation between SNP and traits, the closer Δ(SNPindex) is to 1 [34][35][36]. QTL sequencing was performed by Wuhan Kinosec Technology Co., Ltd.

Simple sequence repeat (SSR) marker development
New SSR markers were developed based on the available genome sequence information of Damor (http://www.genoscope.cns.fr/brassicanapus/) in the candidate region for cqSD-A9a identified by QTL-mapping. The SSR site was obtained through SSRhunter1.3 software, and SSR primers were designed using Primer5.0 software. The primer length ranged from 17-25 bp with an annealing temperature of 50-62˚C, and the length of the amplified fragment was between 100-300 bp (synthesized by Shanghai Sangon Biotech. Co., Ltd.). Polymorphisms were screened in pairwise comparisons using the designed primers. The pairs compared included No. 935 versus No. 3641 and the HSD pool versus the LSD pool. The DH population was used to further narrow the physical interval via the identified polymorphic SSR marker.

RNA sequencing (RNA-seq) analysis
RNA sequencing was conducted for bud (3-5 days before flowering), leaf (three-true-leaf stage), silique (3-5 days after flowering with the same size) and two descendant bud pooled RNA samples. The two descendant pooled samples were prepared by mixing the RNA of 20 high SD DH lines (high SD pool) and 20 low SD DH lines (low SD pool) from the DH population, and all samples were repeated 3 times. The cDNA library was obtained by PCR enrichment after purification by AM Pure XP beads was end-repaired, and A-tailed and sequenced adapters were connected. Clean high-quality data were obtained by filtering and evaluation, were subjected to paired-end 150 bp (PE150) sequencing on the Illumina HiSeq platform. String Tie 2.1.4 (parameters: default) software was used to assemble transcripts, and the merge function of String Tie was applied to the gtf files assembled by each sample to assemble the transcripts obtained from all the samples. These gtf files were merged together using Gff compare 0.12.1 (parameter: -RCK) to compare the merged transcripts with the known transcripts of the genome, discover new transcripts and new genes, and compare existing annotations to obtain supplementary information. We determined the Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways where the differentially expressed genes (DEGs) were significantly enriched relative to all the annotated genes by cluster Profiler 3.14.3 software. RNA-seq analysis of the raw data was performed by staff at Wuhan Bena Technology Service Co., Ltd.B.
All the predicted candidate genes and annotation information within the target interval were obtained from the B. napus Darmor-bzh genome database (http://www.genoscope.cns.fr/ brassicanapus/data/). The Arabidopsis thaliana homologous genes were identified by BLAST analysis of the content within The Arabidopsis Information Resource (TAIR; http://www. arabidopsis.org/), and only the homologous genes that matched the best alignment were screened. The DEGs were summarized and classified to draw a Venn Diagram through-Venny2.1.0 software and subjected to GO functional enrichment analysis to determine the main DEGs enriched in biological processes, cellular components and molecular functions in each tissue. By combining the findings of a KEGG pathway enrichment analysis with the genes identified in the QTL mapping interval, we identified the genes related to SD.

Genetic linkage map construction
A high-density genetic linkage map was constructed using 213 DH lines, which contained 2,102 bins and 1,098,259 single-nucleotide polymorphisms (SNPs) distributed on19 linkage groups, and the total length was 2115.39 cM. The genetic distance of the linkage group ranged from 42.495 cM to 188.675 cM, the average distances between adjacent markers and between adjacent bins were 0.003 cM and 1.058 cM, respectively. The maximum spacing between adjacent markers was 5.389 cM (Table 1 and S1 Table).

Phenotype and correlation analysis
No. 3641, a line with excellent agronomic traits and high SD, was used as the male parent, and No. 935, a line with low SD, was used as the female parent. Furthermore, the SD of No.3641 was significantly higher than that of No. 935 (Fig 1). The differences between the two parents

PLOS ONE
of the DH population for SD were analyzed using t tests and ANOVA. As a result, No.3641 exhibited a higher SD than No.935 in four environments, and the DH lines showed great variation in SD in multiple environments and exhibited high heritability ( Table 2 and S2 Table). Furthermore, the SD performances of the DH lines showed continuous variation and a normal distribution, suggesting that SD is a quantitative trait under polygenic control suitable for QTL analysis (Fig 2). The phenotypic measurements of SD in the four environments of 2019HZ, 2019XN, 2020HZ, and 2020XN were analyzed (Table 3), Correlation analysis revealed that the associations between SD measures across environments were highly significant. This findings indicates that the trait is stable in each environment and less affected by the environment, which will be further verified in the epistasis analysis.
In addition, the phenotypic correlation analysis between SD and other yield-related traits (ESP, SPS and TSW) was performed using SPSS software. The results revealed that (Table 4)

QTL analysis of SD
A total of twenty-eight QTLs for SD were detected in all environments by WinQTLcart 2.5, which were located on chromosomes A02, A04, A05, A09, C02, C03, C06, and C09. The PVE ranged from 4.38% to 13.24%. Eight QTLs for SD located on chromosome A09 (S1 Fig) could explain 5.89%-13.24% of the phenotypic variation, and qSD-A9-1a and qSD-A9-3a were the largest, reaching 13.24% and 12.09%, respectively. In addition, five QTLs located on chromosome C06 could explain 5.61%-11.07% of the phenotypic variation, And 4 QTLs located on chromosome A05 contributed to phenotypic variation at low rates, ranging from 5.96% to 7.92% (Table 5). In summary, SDs were distributed on multiple chromosomes, and the maximum phenotypic variation value was only 13.24%, suggesting that the SD of B. napus is jointly controlled by multiple minor genes. In this study, a total of 17 consistent QTLs for SD were integrated from twenty-eight QTLs by meta-analysis, of which cqSD-A9a was detected in all four environments and located within 7.80-10.43 Mb on chromosome A09, explaining 10.68% of the phenotypic variation, and another consistent QTL, cqSD-C6b, was stably detected in three environments and located on chromosome C06 within the 10-20.41 Mb region, explaining 10.19% of the phenotypic variation. The remaining 3 consistent QTLs were all integrated by less than three environments and located on chromosomes A05, C06, and C09, and all explained less than 10% of the phenotypic

PLOS ONE
variation (Tables 5 and 6). The consistent QTL cqSD-A9a on chromosome A09 showed a positive additive effect, suggesting that the female parent 'No.935' contributed favorable alleles, and the synergistic effect on SD reached 0.13. The contribution rate was the highest and reached 10.68%. Therefore, this locus is considered as major QTL (Table 6).

PLOS ONE
Four pairs of QTL interactions from different chromosomes were detected by QTL epistasis analysis in the DH population: qSD-A02×qSD-C03, qSD-A05×qSD-C07, qSD-C06×qSD-C09, and qSD-A02 ×qSD-A03 (S3 Table). qSD-A02×qSD-C03 showed an (additive × additive) × environment interaction effect and was located on chromosome A02. qSD-A05×qSD-C07 and qSD-C06×qSD-C09 had additive × additive interaction effects; the effect sizes were 1.11% and 1.13%, respectively. qSD-A02×qSD-A03 had no additive effect site, but there was an epistatic effect; the epistatic effect of this interaction was the largest, explaining 3.99% of the phenotypic variation in SD (Fig 3). The QTL × environment interaction had a negative effect, and the effect value was close to 0 in all environments, indicating that there was a certain epistatic interaction for the trait, but it was little affected by the environment.

QTL-seq analysis of SD
To identify the QTL region contributing to SD, QTL sequencing (QTL-seq) was performed, and the Δ (SNP index) was used to measure the allele frequency difference between pools (the B. napus reference genome version used for the alignment is that of Darmor 4.1, available at http://www.genoscope.cns.fr/brassicanapus/data/). The DISTANCE method was used to fit the Δ (SNP index). Then according to the association threshold, the interval above the threshold was considered the candidate region of QTLs for SD. The distribution of the SNP index and Δ (SNP index) of the two mixed pools is shown in S2 Fig. Using the SNP sites with different genotypes between the two pools, we quantified the depth of each base in different pools, and the ED value of each site was calculated. Then, the DISTANCE method was used to fit the ED value, and the distribution of the associated values was subsequently determined, which is shown in Fig 4. According to the distribution of Δ (SNP index) and the association values on the linkage group, a candidate interval of QTL for SD (6.04 Mb to 11.21 Mb) was obtained on chromosome A09. Interestingly, physical interval of the consistent QTL (cqSD-A9a) locus on chromosome A09 detected by QTL mapping ranged from 7,823,737 bp to 10,428,613 bp, which was within the candidate interval of QTLs for SD identified by QTL-seq. It was shown that the consistent QTL for SD (cqSD-A9a) located on chromosome A09 was reliable. In addition, there were a total of 20,764 SNPs and 4,754 insertion-deletions (In Dels) in the candidate interval. The annotation results of these variants are summarized in S4 and S5 Tables.

Development of SSR markers targeting the candidate interval
To confirm the candidate interval of QTLs for SD identified by QTL mapping together with QTL-seq, 18 polymorphic SSR markers (S6 Table) from the candidate region were developed and evaluated in the DH population. As a result, the candidate interval was ultimately narrowed to a 1.86 Mb region (Fig 5).

Prediction of candidate genes related to SD
To help identify the candidate gene and elucidate the molecular basis underlying SD in B. napus, the DEGs were analyzed, and 463 DEGs (308 upregulated, 155 downregulated) were identified between LMPBud and HMPBud, and 1,562 DEGs (515 upregulated, 1047 downregulated) were identified between No.3641Bud and No.935Bud. Similarly, 1731 DEGs (577 upregulated, 1,154 downregulated) were identified in No.3641Leaf and No.935Leaf, and 521 DEGs (48 upregulated, 473 downregulated) were identified in No.3641Pod and No.935Pod. The number of downregulated genes in the different tissues was greater than that of upregulated genes (S3 Fig). Mapping software was used to analyze the intersection and overlap of DEGs, and 485 genes were differentially expressed in different tissues (S4 Fig). A total of 535 DEGs were enriched in 60 different GO terms (Fig 6). In terms of biological processes, the DEGs were mainly enriched in translation, protein ubiquitination, methylation, DNA repair, intracellular protein transport, and transcriptional regulation. Among cellular components, the DEGs were mainly enriched in membrane components, the nucleus, the cytoplasm, the cytosol, ribosomes and chloroplast components. In terms of molecular functions, the DEGs were mainly annotated to ATP binding, DNA binding, RNA binding, nucleic acid binding, oxidoreductase activity, etc. For DEGs identified by KEGG pathway enrichment analysis (S7 Table), several common metabolic pathways were identified in various tissues, including pathways involving selenium compound metabolism, fatty acid biosynthesis, fatty acid metabolism, biotin metabolism, and lysine degradation. We focused on the GO enrichment results of the flower buds and siliques of the parents (S5 and S6 Figs) and found that, in terms of biological processes, the DEGs in these two organs were enriched in translation and DNA repair; cellular components, membrane components, the nucleus, and cytoplasmic components; and molecular function, ATP binding, DNA binding, and metal ion binding.GO enrichment and KEGG analysis showed that DEGs in rapeseed with significant differences in SD are mainly related to protein transport, transcription regulation, DNA binding, etc. The above biological processes play an important role in the formation of seeds. The differential expression of such genes is likely to be an important reason for the differences in SD.
In this study, according to the transcriptome sequencing results, a total of 13 DEGs were screened in the candidate interval identified by QTL mapping (S8 Table). Three of 13 candidate genes were screened out by functional comparison with homologous genes in Arabidopsis. BnaA09g14070D is a gene encodes a callose synthase and was expressed in No. 3641 but not in No. 935 in all tissues. Furthermore, the expression of BnaA09g14070D was also extremely significant between the high SD and low SD pools of the DH population (Fig 7A). Bna09g14800D encodes a membrane component with functions such as intracellular protein transport and vesicle docking, which is homologous to a synaptic protein in Arabidopsis

PLOS ONE
thaliana. The expression difference of Bna09g14800D was extremely significant between the siliques of the two parents ( Fig 7B). BnaA09g18250D functions in DNA binding, transcriptional regulation, and sequence-specific DNA binding and is homologous to the homeoboxleucine zipper genes involved in the metabolic pathway of auxin stimulation in Arabidopsis thaliana. For the expression of BnaA09g18250D, there were extremely significant differences between the buds, leaves and siliques of the two parents as well as the buds of the two extreme pools (extremely high SD and low SD) (Fig 7C).

Discussion
In our study, we evaluated yield-related agronomic traits of No. 935 Table). However, No.935 and No.3641 have no obvious differences in the TSW and SYP. Further research showed that SD had a significant negative correlation with SL and TSW, and a significant positive correlation with SPS and YP via the correlation analysis (Table 4). Some researchers have shown that SL and TSW are important factors affecting the yield of rapeseed. Compared with short siliques, long siliques will produce more and larger seeds. Zhang et al (2011) observed, a positive correlation between SPS and SL, and a significant negative correlation between seed weight and SPS by evaluating Analysis of the differential expression of candidate genes. A represents the relative differential expression of the candidate gene BnaA09g14070D in different tissues, B represents the relative differential expression of the candidate gene Bna09g14800D in different tissues, C represents the relative differential expression of the candidate gene BnaA09g18250D in different tissues. LMPBud: bud mixed pool with extremely low SD; HMPBud: bud mixed pool with extremely high SD, P1: No.935, P2: No. 3641. the silique traits of 140 DH lines and their corresponding parents [10]. Wu et al (2012), showed that appropriately increasing the seed density can increase the yield theoretically with the number of effective siliques, seed size and silique length are constant [37]. In general, the SD, SL, and TSW are mutually restrictive, but they all have a positive effect on increasing yield, as increasing SD can make up for the negative effects of low ESP and TSW traits on yield, Furthermore, increasing SD beneficial for increasing the yield while ensuring a certain seed size and having a suitable silique length and shape.
At present, there are few studies on the SD of spring B. napus. Only a few researchers have performed preliminary research on this trait, and several related QTLs and genes have been found. Ren et al., investigated the phenotype of SD and its related traits in natural populations and found that several traits were normally distributed. Two SNPs were found on chromosomes A07 and A10 via GWAS [4]. The SD-related QTLs discovered are distributed on chromosomes C04, C06 and C09 [14]. The SD QTLs detected were previously located on chromosomes A09 and C06 by Wang et al. [38]. Li [5] detected two loci that control the number of grains per unit length; these loci are located in linkage groups N6 and N19. Deng et al. [26] detected 13 target QTLs related to SD on A3, A7, A9 and C3, and 3 important QTLs were located at 26.58-29.19 Mb on chromosome A09. In our study, a high-density genetic linkage map was constructed by resequencing. A total of 28 QTLs for SD were detected by Win QTL cart 2.5, which were located on chromosomes A02, A04, A05, A09, C02, C03, C06, and C09. The consistent QTL cqSD-A9a on chromosome A09 showed a positive additive effect with a 10.68% contribution rate and was considered the major QTL. The candidate interval of QTLs for SD was ultimately narrowed to a 1.86 Mb region on A09 by QTL mapping together with QTL-seq and SSR marker encryption. On the basis of the above conclusions, it can be seen that the QTLs that control SD are distributed on multiple chromosomes. Although the major QTL was located on chromosome A09, it is not in the range of previous studies, and in this study, it is likely to be a new site. We can further identify genes within this candidate interval.
To date, research on SD-related genes is still limited, and only a few related functional genes have been screened by individual researchers. These functional genes are related to the regulation of cytokinins, the formation process of seeds and the growth and development of lateral organs [4]. In this study, according to the transcriptome sequencing results, 13 DEGs were screened in the candidate interval identified by QTL mapping, and 3 of 13 candidate genes were screened out by functional comparison with homologous genes in Arabidopsis. BnaA09g14070D (AT5G13000) encodes a callose synthase-like enzyme that plays an important role in developmental and stress responses [39]. Callose is a β-1,3-bonded glucan that protects mother cells from external stimuli and harmful substances during gametophyte development, and is of great significance to the development of male and female gametes. BnaA09g14800D (AT3G09740), also known as ATSYP71 (plant synapsin 71), mediates membrane fusion in Arabidopsis cytokinesis [40]. Synapsin belongs to the Qa-SNARE (soluble N-ethylmaleimidesensitive factor linking complex) family and participates in various biological processes [41,42]. Plant Qa-SNARE proteins also membrane material transport, affecting cell division and plant growth differentiation [43]. BnaA09g18250D (AT5G47370) is also referred to as the homeobox leucine zipper 4 (HB-4)/HD-ZIP protein. Basic leucine zipper (bZIP) transcription factors control important developmental and physiological processes in plants [44]. The Znregulated transporter iron-regulated transporter (IRT)-like protein (ZIP) family members are involved in Zn transport and cellular Zn homeostasis throughout the domains of life. IRT3, ZIP4, ZIP6, and ZIP9 function redundantly in maintaining Zn homeostasis and seed development in A. thaliana [45]. The functions of these candidate genes are all related to the regulation of cytokinins, the formation process of seeds and the growth and development of lateral organs. Whether these genes control the SD of rapeseed seeds remains to be verified.
In summary, these findings are important because they provide information regarding seed density per silique (SD) in B. napus. They also provide a basis for future breeding and gene cloning.

Conclusion
In the present study, a genetic linkage map was constructed containing 2,102 bins and 1,098,259 single-nucleotide polymorphisms (SNPs) distributed among 19 linkage groups via a DH population of 213 lines derived from the F 1 cross between two varieties with significant differences in SD. A major QTL cqSD-A9a was mapped within 7.80 Mb-10.43 Mb on chromosome A09 controlling the SD trait by QTL-mapping and QTL-Seq analysis. Three differentially expressed genes, i.e., BnaA09g14070D, BnaA09g14800D, and, BnaA09g18250D may be associated with SD and were obtained in the candidate interval by RNA-Seq analysis.